//
// Created by tview on 20-12-22.
//


struct Complex  {
	double r, i ;
	Complex ( )  {	}
	Complex ( double r, double i ) : r ( r ), i ( i )  {	}
	inline void real ( const double& x )  {  r = x ;  }
	inline double real ( )  {  return r ;  }
	inline Complex operator + ( const Complex& rhs )  const  {
		return Complex ( r + rhs.r, i + rhs.i ) ;
	}

	inline Complex operator - ( const Complex& rhs )  const  {
		return Complex ( r - rhs.r, i - rhs.i ) ;
	}

	inline Complex operator * ( const Complex& rhs )  const  {
		return Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
	}

	inline void operator /= ( const double& x )   {
		r /= x, i /= x ;
	}

	inline void operator *= ( const Complex& rhs )   {
		*this = Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
	}

	inline void operator += ( const Complex& rhs )   {
		r += rhs.r, i += rhs.i ;
	}

	inline Complex conj ( )  {
		return Complex ( r, -i ) ;
	}
} ;

bool inverse = false ;

inline Complex omega ( const int& n, const int& k )  {
    if ( ! inverse ) return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ) ;
    return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ).conj ( ) ;
}

//木的优化的版本
inline void fft ( Complex *a, const int& n )  {
    if ( n == 1 ) return ;

    static Complex buf [N] ;

    const int m = n >> 1 ;

    for ( int i = 0 ; i < m ; ++ i )  {
        buf [i] = a [i << 1] ;
        buf [i + m] = a [i << 1 | 1] ;
    }

    memcpy ( a, buf, sizeof ( Complex ) * ( n + 1 ) ) ;

    Complex *a1 = a, *a2 = a + m;
    fft ( a1, m ) ;
    fft ( a2, m ) ;

    //快速傅立叶变换根本原因是分治，也就是相邻两个系数是可以通过共同的计算数值。不要纠结着绕的算法，大概记忆一下。
    for ( int i = 0 ; i < m ; ++ i )  {
        Complex t = omega ( n, i ) ;
        buf [i] = a1 [i] + t * a2 [i] ;
        buf [i + m] = a1 [i] - t * a2 [i] ;
    }

    memcpy ( a, buf, sizeof ( Complex ) * ( n + 1 ) ) ;
}
struct FastFourierTransform  {
    Complex omega [N], omegaInverse [N] ;

    void init ( const int& n )  {
        for ( int i = 0 ; i < n ; ++ i )  {
            omega [i] = Complex ( cos ( 2 * PI / n * i), sin ( 2 * PI / n * i ) ) ;
            omegaInverse [i] = omega [i].conj ( ) ;
        }
    }

    void transform ( Complex *a, const int& n, const Complex* omega ) {
        //换位，本来是通过递归，但是这里通常使用循环实现。
        for ( int i = 0, j = 0 ; i < n ; ++ i )  {
		    if ( i > j )
		        std :: swap ( a [i], a [j] ) ;
		    //ｌ＝ｎ/2          ｊ＝ｌ＝》　ｊ＝ｊ取最位                 l = l/2
		    for( int l = n >> 1 ; ( j ^= l ) < l ; l >>= 1 ) ;
	    }

        for ( int l = 2 ; l <= n ; l <<= 1 )  {
            int m = l / 2;
            for ( Complex *p = a ; p != a + n ; p += l )  {
                for ( int i = 0 ; i < m ; ++ i )  {
                    Complex t = omega [n / l * i] * p [m + i] ;
                    p [m + i] = p [i] - t ;
                    p [i] += t ;
                }
            }
        }
    }

    void dft ( Complex *a, const int& n )  {
        transform ( a, n, omega ) ;
    }

    void idft ( Complex *a, const int& n )  {
        transform ( a, n, omegaInverse ) ;
        for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ;
    }
} fft ;




/*
*做FFT
*len必须是2^k形式
*on == 1时是DFT，on == -1时是IDFT
*/
void fft(Complex y[],int len,int on){
    change(y,len);
    for(int h = 2;h <= len;h<<=1){
        Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0;j < len;j += h){
            Complex w(1,0);
            for(int k = j;k < j + h/2;k++){
                Complex u = y[k];
                Complex t = w*y[k + h/2];
                y[k] = u + t;
                y[k + h/2] = u - t;
                w = w*wn;
            }
        }
    }
    if(on == -1){
        for(int i = 0;i < len;i++){
            y[i].x /= len;
        }
    }
}


/*
*做FFT
*len必须是2^k形式
*on == 1时是DFT，on == -1时是IDFT
*/
void fft(Complex y[],int len,int on){
    change(y,len);
    for(int h = 2;h <= len;h<<=1){//模拟合并过程
        Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));//计算当前单位复根
        for(int j = 0;j < len;j += h){
            Complex w(1,0);//计算当前单位复根
            for(int k = j;k < j + h/2;k++){
                Complex u = y[k];
                Complex t = w*y[k + h/2];
                y[k] = u + t;//这就是吧两部分分治的结果加起来
                y[k + h/2] = u - t;
                //后半个“step”中的ω一定和“前半个”中的成相反数
                //“红圈”上的点转一整圈“转回来”，转半圈正好转成相反数
                //一个数相反数的平方与这个数自身的平方相等
                w = w*wn;
            }
        }
    }
    if(on == -1){
        for(int i = 0;i < len;i++){
            y[i].x /= len;
        }
    }
}